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ABSTRACT 


The  accuracy  of  the  SPACETRACK  General  Perturbations  progi'am  over 
short  periods  is  evaluated  for  a  number  of  cases.  They  show  that  oscillatory 
terms  resulting  from  drag  perturbations  contribute  heavily  to  errors  at  low 
altitudes,  and  that  these  terms  must  be  eliminated  if  the  first  order  theory 
is  to  be  used  for  high  accuracy  in  these  circumstances.  It  appears  that  the 
evaluation  technique  employed  would  be  useful  in  addressing  a  number  of 
other  problems;  several  promising  applications  are  discussed. 
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GLOSSARY 


A  cross-sectional  area 

a  semimajor  axis 

drag  coefficient 
e  eccentricity 

H  atmospheric  scale  height 

h  perigee  in  distance  above  earth's  surface 

L  mean  longitude 

m  mass  of  satellite 

—  n  —  rate  of  change  of  mean  motion 

2  2 

P  period 

q  perigee  in  distance  from  center  of  earth 

t  time  since  epoch 

v  velocity  of  satellite  with  respect  to  air  mass 

/?  ballistic  coefficient 

p  atmospheric  density 

atmospheric  density  at  perigee 
X  dimensionless  drag  parameter  =  p^  0  q 


SECTION  I 


INTRODUCTION 

In  the  SPACETRACK  System,  it  is  customary  to  use  a  special  pertur- 

bations  propram  when  a  high  depree  of  accuracy  is  desired  over  a  short 

interval  of  time.  The  purpose  of  this  investipation  is  to  determine  whether 

*  * 

and  under  what  circumstances  a  first  order  peneral  perturbations  propram 
may  be  used  instead;  and  if  possible,  determine  whether  any  modifications 
to  it  are  feasible.  The  interval  of  time  considered  in  this  study  is  of  the 
order  of  1-1/2  days,  and  thus  interest  is  in  accuracy  of  the  order  of  1  km. 
Most  of  the  study  deals  with  low  altitudes  since  drap  is  the  perturbative 
force. 

There  already  was  in  existence  before  the  start  of  this  investipation  a 
propram  called  DCMOD64,  written  by  the  Aeronutronic  Division  of  the  Philco 
Corporation,  well  suited  for  the  purpose  of  makinp  comparisons.  In  fact, 
without  this  propram's  prior  existence,  it  would  have  been  impossible  to 
finish  this  study  in  a  reasonable  time.  Other  investipators  have  already  used 
this  propram  for  similar  studies.^  In  this  study,  we  are  concerned  with  a 
shorter  period  for  the  approximation  interval  and  hipher  accuracies. 


A  program  in  which  the  effects  of  the  perturbing  forces  are  numerically 
integrated. 

A  program  which  employs  an  analytic  theory  of  the  effects  of  perturbing 
forces. 
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SECTION  II 


DESCRIPTION  OF  COMPUTER  PROGRAM 

[2] 

The  DCMOD64  program  includes  an  orbital  element  correction  routine, 

and  several  ephemeris  computation  subroutines.  We  have  used  a  special 

option  of  this  program  for  our  study.  In  this  option,  the  special  perturbations 

subroutine  generates  an  ephemeris  from  an  initial  set  of  orbital  elements. 

This  is  then  converted  into  400  equally  spaced  (in  time)  observations  from  a 

hypothetical  radar  with  spherical  coverage.  Then  the  general  perturbations 

subroutine  is  used  in  the  element  correction  routine  to  fit  these  observations. 

The  final  output  is  a  table  of  discrepancies  between  the  special  and  general 

perturbations  ephemerides.  The  first  order  general  perturbations  theory 

is  equivalent  to  Lyddane's  modification  of  the  Brouwer  theory;  it  should  achieve 

6  3 

precisions  on  the  order  of  1  part  in  10  ,  or  1  microradian,  over  10  radians 
of  satellite  motion,  insofar  as  perturbations  due  to  the  earth's  potential  are 
concerned.  Only  the  3  largest  zonal  harmonics,  J  -  J  ,  are  included  in  the 
formulation.  The  subroutine  also  includes  a  formulation  for  the  perturbative 
effects  of  solar  radiation  pressure,  which  were  of  little  significance  in  this 
study.  Other  perturbations,  of  which  air  drag  is  the  most  significant,  are 
accommodated  by  two  empirical  terms  in  the  mean  anomaly  equation,  n/2 
and  n/6  ,  so  that  a  correction  to  the  mean  anomaly  is  given  in  the  form 

1.2  1  ..  3 

6M  =  —  n  t  +  —  nt  , 

2  6 


where  t  is  time  since  epoch. 

Related  corrections  to  the  semimajor  axis  and  the  eccentricity  are  derived 
from  these  parameters  under  the  assumption  of  constant  perigee  height.  The 
special  perturbation  subroutine  that  generates  the  ephemeris  takes  into  account 


2 


eight  zonal  harmonics,  four  tesseral  harmonics,  atmospheric  drag,  solar 

radiation  pressure,  and  lunar  and  solar  gravitational  perturbations.  It  is 

possible  to  omit  any  or  all  of  these  in  any  rim.  The  atmosphere  model  takes 

into  account  the  diurnal  bulge  as  well  as  solar  activity,  but  these  are  not 

effective  at  the  altitudes  chosen  for  this  study.  The  ballistic  coefficient,  /?  , 

is  assumed  constant.  For  the  puspose  of  this  study,  it  was  necessary  to 

modify  the  DCMOD  program  with  octal  corrections  to  make  the  time  period 
* 

36  hours. 


The  authors  are  indebted  to  the  cooperation  and  programming  assistance  of 
J.  Kuhlman  of  Aeronutronic  for  obtaining  data  successfully. 
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SECTION  HI 


DESCRIPTION  OF  CASES 


Several  different  classes  of  satellites  were  simulated.  All  were  specified 
by  their  initial  osculating  elements.  All  had  an  initial  inclination  of  49  degrees. 
Initial  perigees  used  were: 


h  (km)  h  (km) 


160 

87.  7 

200 

109.6 

250 

137.  0 

300 

164.4 

These  were  used  in  combination  with  initial  eccentricities  of: 

0.  0 
0.  001 
0.  01 
0.  1 
0.2 
0.  3 


Various  combinations  of  the  perigees  and  eccentricities  were  used  in  the 
program  under  three  different  classes  of  perturbations  in  the  special  pertur¬ 
bations  program: 

(a)  all  perturbations; 

(b)  only  second,  third,  and  fourth  zonal  harmonics  of  earth's  potential; 
and 

(c)  only  atmospheric  drag  and  second,  third,  and  fourth  zonal  harmonics 
of  earth's  potential. 
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It  was  originally  intended  that  all  satellites  have  the  same  ballistic 

2  . 

coefficient,  p  .  The  value  0.02  m  /kg  was  tried,  since  this  appears  to  be 

[3] 

a  high  average,1  and  results  using  this  value  should  be  conservative. 
However,  for  some  low-altitude  satellites  with  small  eccentricities,  the 
special  perturbations  subroutine  would  not  run  36  hours  with  this  value,  pre¬ 
sumably  because  of  decay,  and  so  the  smaller  values  listed  in  the  tables  were 
used  in  these  cases. 
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SECTION  IV 


SUMMARY  OF  RESULTS 

Thirty-six  different  cases  were  run;  the  results  are  summarized  in  Tables 

I,  II,  and  III.  Each  of  these  tables  presents  the  distinguishing  initial  conditions 

* 

followed  by  columns  containing  five  mean  parameters  determined  by  AGP; 
viz,  h  (computed  perigee),  e  (eccentricity),  a  (semimajor  axis),  P  (period), 
and  n/2  (rate  of  change  of  mean  motion).  These  are  followed  by  the  errors 
in  36  hours  of  simulated  time.  The  RMS  error  is  the  root  mean  square  of  all 
400  error  vector  magnitudes  in  the  36-hour  period.  The  maximum  error  is 
the  maximum  magnitude  of  the  400  vector  errors  in  the  same  time  period. 

Table  I  contains  the  cases  wherein  all  perturbations  are  included.  Table 
II  presents  the  cases  wherein  the  only  perturbations  are  the  2nd,  3rd,  and 
4th  zonal  harmonic  terms  of  the  earth's  potential.  Table  IIT,  which  presents 
the  cases  in  which  the  only  perturbations  are  atmospheric  drag  and  the  oblate 
earth  terms,  contains  two  additional  columns,  the  first  of  which  is  adjusted 
RMS  errors.  These  are  the  RMS  values  adjusted  by  multiplying  the  factor 
0.  02//3  .  This  is  done  to  facilitate  comparison  of  results  for  different  cases 
since,  presumably,  errors  are  roughly  proportional  to  the  ballistic  coefficient. 
The  figures  in  the  other  additional  column  are  normalized  RMS  errors.  These 
are  nondimensional  quantities  consisting  of  the  RMS  error  value  divided  by 
the  quantity  (ax/e),  where  x-  Pq/3  q  is  a  nondimensional  parameter  which 
gives  the  order  of  magnitude  of  the  errors  due  to  atmospheric  drag.  A  fuller 
discussion  of  these  quantities  is  presented  later. 


* 

AGP  is  the  acronym  of  the  general  perturbations  routine  designed  for  non- 
equatorial  cases. 
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Errors  in  Approximating  Ephemris  by  General  Perturbations  - 
(Only  atmospheric  and  2nd,  3rd,  and  4th  zonal  harmonic  perturbations  included) 
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A  point  of  interest  that  does  not  appear  in  these  tables  is  the  nature  of 
the  error  as  a  function  of  time,  or  alternatively,  as  a  function  of  true  agrument 
of  latitude.  Generally  speaking,  in  the  cases  in  which  only  atmospheric  and 
zonal  harmonic  perturbations  are  considered,  the  out-of-plane  component  is 
an  order  of  magnitude  smaller  than  the  other  two  components.  Both  the 
along-track  and  radial  components  have  a  decidedly  oscillatory  behavior, 
with  maximum  amplitude  at  the  ends  of  the  simulation  interval  and  a  phase 
change  near  the  center  of  the  interval.  Figure  1  shows  the  three  components 
of  the  error  as  a  function  of  observation  number,  which  is  proportional  to 
time,  for  job  213,  which  may  be  regarded  as  archetypical  of  the  runs  in 
Table  III. 
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250  KM 
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SECTION  V 


DISCUSSION  OF  RESULTS 


GENERAL 

The  runs  in  Table  I  were  made  primarily  to  determine  regions  wherein 
the  first-order  general  perturbations  theory  is  adequate  by  comparison  with 
the  best  available  model  of  the  real  world.  They  clearly  show  that  for  satel¬ 
lites  with  a  160-km  perigee,  the  theory  will  give  errors  much  larger  than 
1  km  for  eccentricities  less  than  0.01.  A  comparison  with  the  cases  in 
Table  III  at  this  same  perigee  shows  a  very  high  correlation.  This  strongly 
suggests  that,  at  this  perigee,  the  oblate  earth  and  atmospheric  perturbations 
are  the  only  ones  which  are  significant  for  short  periods.  This  was  the  reason 
for  confining  the  rest  of  the  study  to  an  examination  of  drag  effects.  A  com¬ 
parison  of  the  data  at  the  300-km  perigee  shows  a  much  poorer  correlation  — 
other  perturbations  are  much  more  significant  in  proportion  —  but  the  others 
would  still  be  small  enough  to  be  acceptable  if  the  drag  errors  could  be 
eliminated.  In  general,  the  data  in  Table  I  suggest  that  the  existing  first- 
order  theory  is  adequate  for  any  altitude  above  350  km. 

The  runs  in  Table  II,  with  only  J  -  J  zonal  harmonic  perturbations, 

U4  fr 

were  intended  to  provide  a  comparison  for  the  runs  in  Table  III  which  has 
atmospheric  drag  as  well  as  the  J  -  J  zonal  harmonic  perturbations. 

Ud  ft 

Since  the  first-order  theory  accounts  for  these  zonal  harmonic  perturbations, 
the  errors  should  be  on  the  order  of  6  to  18  meters  for  all  Table  II  cases. 

The  reason  they  are  not  zero  is  not  entirely  clear. 

Part  of  the  difference  arises  from  the  fact  that  the  J  -  J.  values 

2  4 

stored  in  the  special  and  general  perturbations  subroutines  do  not  agree. 

This  disagreement  arises  because  the  special  perturbations  set,  with 
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12  parameters,  and  the  general  perturbations  set,  with  3  parameters,  have 
been  independently  adjusted  for  a  best  fit  with  the  observed  motion  of  satellites . 
This  discrepancy  was  not  discovered  sufficiently  early  in  the  study  for  cor¬ 
rective  measures  to  be  taken.  As  a  check  on  the  significance  of  the  discrepancy, 
job  115  was  rerun.  The  maximum  error  was  reduced  from  308  meters  to 
18  meters,  while  the  RMS  error  was  reduced  from  68  meters  to  6.  8  meters. 

The  oscillatory  pattern  of  the  errors  remains,  but  there  is  no  apparent 
tendency  for  the  amplitude  of  the  oscillation  to  grow  with  time. 

A  second  phenomens  is  evident  in  Table  II:  the  errors  grow  with  de¬ 
creasing  eccentricity.  It  is  possible  that  this  is  due  to  the  discrepancy  in 
the  J  -  J  terms;  on  the  other  hand,  it  may  reflect  numerical  problems 

Z  A 

in  the  element  correction  process. 

The  convergence  of  the  solutions  appeared  to  be  rather  slow  in  all  cases; 

* 

from  6  to  8  Phase  II  iterations  were  usually  required  for  convergence  to 
a  1 -percent  change  in  the  RMS  of  the  vector  magnitudes. 

DRAG  EFFECTS 

Table  III  presents  the  errors  in  runs  with  drag  perturbations  as  well  as 
the  J  -  J  zonal  harmonics.  Since  the  general  perturbations  subroutines 
account  for  the  harmonics  with  the  accuracies  given  in  Table  II,  the  additional 
errors  in  these  cases  must  be  dxie  to  drag  alone  or  to  cross-coupling  between 
drag  and  zonal  harmonics.  The  magnitude  of  the  acceleration  of  a  satellite 
due  to  drag  is  given  by 


In  Phase  I,  only  the  mean  anomaly  or  ’’time”  equation  is  corrected;  in 
Phase  II,  all  elements  are  corrected. 
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where  p  is  the  density,  v  is  the  velocity  with  respect  to  the  air  mass,  and 
f3  ,  the  ballistic  coefficient,  is  given  by 


13  =  CD  A/m  , 


where  is  the  drag  coefficient,  A  is  the  cross-sectional  area  of  the 
satellite,  and  m  is  its  mass. 


In  general,  as  a  satellite  rotates,  neither  nor  A  remains  constant 


so  that  (3  varies  with  the  orientation.  However,  for  satellites  in  which  the 
ratio  of  the  longest  to  shortest  dimension  is  no  larger  than  two  or  three,  the 
product  does  not  vary  very  much,  and,  for  a  given  satellite,  /?  varies  per- 


one  satellite  to  another  because  of  differences  of  mass.  The  value  of  0.  02 
2 

m  /kg  used  in  the  simulation  is  conservative  in  that  few  satellites  could  be 
expected  to  have  a  larger  value. 

Air  density,  p  ,  decreases  quite  rapidly  with  altitude.  For  altitudes 
between  160  and  300  km,  the  scale  height,  H  ,  given  by 


[5] 


is  of  the  order  of  25  to  50  km. 


Since  perigee  does  not  change  very  rapidly  for  satellites  above  200  km, 
a  useful  dimensionless  parameter  that  given  the  scale  of  perturbations  is 


X  =  PQ  13  q  , 


where  p  is  density  at  perigee,  q  . 


o 


Pages  12  to  17. 
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The  last  column  of  Table  III  shows  a  definite  correlation  between  RMS 
errors  and  ax/e.  The  reason  for  the  inverse  variation  of  the  RMS  errors 
with  eccentricity  is  not  known.  Quite  possibly  the  problem  will  prove  to  be 
identical  to  that  in  the  nondrag  cases  of  Table  II;  the  amplitude  and  growth 
of  the  oscillations,  however,  is  substantially  greater. 

[6] 

In  an  analysis  in  terms  of  coordinates,  Geyling  obtained  terms  similar 
to  this.  His  factor  of  proportionality  does  not  appear  to  vary  inversely  as 
the  eccentricity;  perhaps  this  is  because  his  atmospheric  density  model  is 
independent  of  altitude.  Analysis  along  these  lines,  using  a  more  realistic 
density  model,  might  show  results  similar  to  those  obtained  here,  although, 
at  best,  preliminary  analysis  has  shown  a  variation  inversely  as  the  square 
root  of  the  eccentricity. 
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SECTION  VI 


CONCLUSIONS 

The  techniques  employed  for  this  study  have  considerable  potential  for 
studying  the  performance  of  the  DCMOD64  element  correction  routines.  The 
rather  slow  convergence  and  possible  poor  performance  for  low  eccentricity 
deserve  further  study.  Other  areas  of  interest  include  the  dependence  of 
quality  of  fit  on  amount  of  data  and  length  of  arc.  By  using  an  identical 
ephemeris  subroutine  for  data  generation  and  fitting,  it  is  possible  to  check 
numerical  and  partial  derivative  problems.  By  using  special  and  general 
perturbations  routines  with  identical  harmonics,  it  is  possible  to  crosscheck 
the  mathematical  formulations,  to  determine  the  intervals  over  which  the 
routines  remain  valid  and  to  determine  what  length  of  fit  is  necessary  to 
prevent  second  order  terms  in  the  semimajor  axis  from  propagating  into 
the  mean  motion.  By  using  the  full  set  of  harmonics  in  the  special  pertur¬ 
bations  subroutine  and  the  J  -  J  set  in  the  general  perturbations  sub- 

Cl  T 

routine,  the  need  for  additional  general  perturbations  formulations  can  be 
assessed. 

Similar  tests  can  be  made  for  the  tesserals  and  drag  perturbations.  In 
the  case  of  drag  perturbations,  it  may  well  be  the  optimal  procedure  to 
develop  additional  general  perturbations  formulations  based  on  an  empirical 
analysis  of  the  periodic  residuals  (for  periodic  terms),  and  on  the  analysis  of 
mean  elements  for  overlapping  arcs  (for  secular  terms).  Among  the  more 
obvious  questions  that  can  be  resolved  by  such  techniques  are  the  extent  of 
gravitational  drag  cross-coupling  in  the  motion  of  node  and  perigee,  and 
the  secular  behavior  of  eccentricity. 

For  greater  consistency  and  to  avoid  errors  in  using  subsets  of  the  full 
potential  model,  it  is  suggested  that  the  DCMOD64  control  logic  for  the 
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special  perturbations  subroutine  be  modified  to  permit  selection  of  optimized 

subsets  of  the  zonal  and  tesseral  harmonics,  rather  than  selective  inclusion 

or  deletion  of  individual  harmonics.  Based  on  reports  in  the  literature  of 

the  techniques  used  to  determine  the  harmonics,  it  is  probably  possible  to 

treat  even  zonals,  odd  zonals,  low-order  tesserals  (n  ^  8),  and  high-order 

tesserals  (n  13)  as  four  independently  optimized  subjects.  Suggested 

options  would  include  the  general  perturbations  subroutine  values  of  JQ  - 

the  current  full  set  of  zonals  (J  -  J  or  J  -  J  )  ,  and  the  full  set  of 

2  9  2  14 

zonal  harmonics  plus  all  available  tesserals. 
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